Wheelset

H. G. Kwatny

Department of Mechanical Engineering & Mechanics
Drexel University

Introduction

<< ProPac`;

?ProPac`Dynamics`*

Wheelset Model

Main body joint

r1 = {3}; q1 = {ϕ, x, y}; p1 = {ωb, vx, vy};
H1 = {{0, 0, 0}, {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 
            0}};

JointLst = {{r1, H1, q1, p1}};

Body data

cm1 = {0, 0, 0}; 
    out1 = {{2, {0, d, 0}}, {3, {0, 0, 0}}, {4, {0, -d, 
            0}} } ;   I1 = {{Jxx, 0, Jxz}, {0, Jyy, 0}, {Jxz, 0, 
          Jzz}};

BodyLst = {{cm1, out1, m1, I1}};
TreeLst = {{{1, 1}}};
q = q1;
p = p1;

Q = {T, F, 0};
\(\({V, X, H, M, Csys, Fsys, psys, qsys} = 
	
        CreateModel[JointLst, BodyLst, TreeLst, g, 0, Q];\)\)

Constraints

Vset = NodeVelocity[3, TreeLst, BodyLst, X, H, q, p];
c1 = Vset[[5]]

p

\(\({Mm, Cm, Fm, Vm, Trans, phat} = 

        Simplify[
          DifferentialConstraints[M, Csys, Fsys, V, {c1}, p, 
            q, {1, 2}]];\)\)

Vm // MatrixForm

q

Simplify[Mm]

f = Simplify[(Cm . phat + Fm)]

Vm . phat

phat

RHS = Simplify[Join[(Vm . phat), -f] /. d -> 1]

Controller & Simulations

Potential Function

U = -Log[1/(Sqrt[x^2 + y^2])];
Plot3D[U, {x, -1, 1}, {y, -1, 1}, 
      AxesLabel -> {x, y, "U"}, PlotPoints -> 25]

EE = -Jacob[U, {x, y}]

f = 
      FullSimplify[EE, x ∈ Reals && y ∈ Reals]

AA = 
      FullSimplify[ArcTan[f[[1]], f[[2]]], 
        x ∈ Reals && y ∈ Reals]

BB = 
      FullSimplify[Jacob[AA, {x, y}], 
        x ∈ Reals && y ∈ Reals]

FullSimplify[BB . f, 
      x ∈ Reals && y ∈ Reals]

s2 = λ (θ - AA) + ω

Reduce[{vx Cos[ϕ] == EE[[1]], 
        vx Sin[ϕ] == EE[[2]]}, {ϕ}, Reals]

Controller

s1 = 
        vx - vbar /. 
          vbar -> (1/
                  10) (Abs[x] + Abs[y])/(1 + Abs[x] + 
                    Abs[y]);
s2 = (ωb + λ (ϕ - 
                    ArcTan[-x, -y])) /. λ -> 
            5;

X = Join[q, phat]

RHScl = (RHS /. {F -> -Sign[s1], 
            T -> -4 Sign[s2]})

eqns = MakeODEs[X, RHScl, t]

Simulations 1

InitialConditions = {x[0] == 0, 
          y[0] == 2/10, ϕ[0] == 0, 
          vx[0] == 0, ωb[0] == 
            0};
Eqns = Join[eqns, InitialConditions];

sols = 
      NDSolve[Eqns, X, {t, 0, 50.0}, AccuracyGoal -> 12, 
        PrecisionGoal -> 12, WorkingPrecision -> 12, 
        MaxSteps -> 50000, MaxStepSize -> 0.01, 
        Method -> {FixedStep, 
            Method -> ExplicitRungeKutta}]

Plot[Evaluate[y[t] /. sols], {t, 0, 50}]

Plot[Evaluate[vx[t] /. sols], {t, 0, 50}, 
      PlotRange -> All]

Plot[
      Evaluate[((ωb[
                    t] + λ (ϕ[t] - 
                        ArcTan[-x[t], -y[
                              t]])) /. λ -> 5) /. 
          sols], {t, 0, 50}, PlotRange -> All]

Plot[Evaluate[ϕ[t] /. sols], {t, 0, 50}, 
      PlotRange -> All]

gr1 = 
        ParametricPlot[Evaluate[{x[t], y[t]} /. sols], {t, 0, 50}, 
          PlotPoints -> 1000, PlotRange -> All, 
          AxesLabel -> {"x", "y"}];

Simulations 2

InitialConditions = {x[0] == 2/10, 
          y[0] == 0, ϕ[0] == 0, 
          vx[0] == 0, ωb[0] == 
            0};
Eqns = Join[eqns, InitialConditions];

Tfinal = 50;
sols = 
      NDSolve[Eqns, X, {t, 0, Tfinal}, AccuracyGoal -> 12, 
        PrecisionGoal -> 12, WorkingPrecision -> 12, 
        MaxSteps -> 50000, MaxStepSize -> 0.01, 
        Method -> {FixedStep, 
            Method -> ExplicitRungeKutta}]

Plot[Evaluate[y[t] /. sols], {t, 0, Tfinal}]

Plot[Evaluate[vx[t] /. sols], {t, 0, Tfinal}, PlotRange -> All, 
      AxesLabel -> {"t", "vx"}]

Plot[
      Evaluate[((ωb[
                    t] + λ (ϕ[t] - 
                        ArcTan[-x[t], -y[
                              t]])) /. λ -> 5) /. 
          sols], {t, 0, Tfinal}, 
      PlotRange -> All . 
            AxesLabel -> {"t", "s2"}]

Plot[Evaluate[ϕ[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

gr2 = 
        ParametricPlot[Evaluate[{x[t], y[t]} /. sols], {t, 0, Tfinal}, 
          PlotPoints -> 1000, PlotRange -> All, 
          AxesLabel -> {"x", "y"}];

Show[gr1, gr2]

Simulations 3

InitialConditions = {x[0] == 2/10, 
          y[0] == -2/10, ϕ[0] == 0, 
          vx[0] == 0, ωb[0] == 
            0};
Eqns = Join[eqns, InitialConditions];

Tfinal = 50;
sols = 
      NDSolve[Eqns, X, {t, 0, Tfinal}, AccuracyGoal -> 12, 
        PrecisionGoal -> 12, WorkingPrecision -> 12, 
        MaxSteps -> 50000, MaxStepSize -> 0.01, 
        Method -> {FixedStep, 
            Method -> ExplicitRungeKutta}]

Plot[Evaluate[y[t] /. sols], {t, 0, Tfinal}]

Plot[Evaluate[vx[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

Plot[
      Evaluate[((ωb[
                    t] + λ (ϕ[t] - 
                        ArcTan[-x[t], -y[
                              t]])) /. λ -> 5) /. 
          sols], {t, 0, Tfinal}, PlotRange -> All]

Plot[Evaluate[ϕ[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

gr3 = 
        ParametricPlot[Evaluate[{x[t], y[t]} /. sols], {t, 0, Tfinal}, 
          PlotPoints -> 1000, PlotRange -> All, 
          AxesLabel -> {"x", "y"}];

Show[gr1, gr2, gr3]

Simulations 4

InitialConditions = {x[0] == -2/10, 
          y[0] == -2/10, ϕ[0] == 0, 
          vx[0] == 0, ωb[0] == 
            0};
Eqns = Join[eqns, InitialConditions];

Tfinal = 50;
sols = 
      NDSolve[Eqns, X, {t, 0, Tfinal}, AccuracyGoal -> 12, 
        PrecisionGoal -> 12, WorkingPrecision -> 12, 
        MaxSteps -> 50000, MaxStepSize -> 0.01, 
        Method -> {FixedStep, 
            Method -> ExplicitRungeKutta}]

Plot[Evaluate[y[t] /. sols], {t, 0, Tfinal}]

Plot[Evaluate[vx[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

Plot[
      Evaluate[((ωb[
                    t] + λ (ϕ[t] - 
                        ArcTan[-x[t], -y[
                              t]])) /. λ -> 5) /. 
          sols], {t, 0, Tfinal}, PlotRange -> All]

Plot[Evaluate[ϕ[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

gr4 = 
        ParametricPlot[Evaluate[{x[t], y[t]} /. sols], {t, 0, Tfinal}, 
          PlotPoints -> 1000, PlotRange -> All, 
          AxesLabel -> {"x", "y"}];

Show[gr1, gr2, gr3, gr4]

Simulations 5

InitialConditions = {x[0] == -2/10, 
          y[0] == 2/10, ϕ[0] == 0, 
          vx[0] == 0, ωb[0] == 
            0};
Eqns = Join[eqns, InitialConditions];

Tfinal = 50;
sols = 
      NDSolve[Eqns, X, {t, 0, Tfinal}, AccuracyGoal -> 12, 
        PrecisionGoal -> 12, WorkingPrecision -> 12, 
        MaxSteps -> 50000, MaxStepSize -> 0.01, 
        Method -> {FixedStep, 
            Method -> ExplicitRungeKutta}]

Plot[Evaluate[y[t] /. sols], {t, 0, Tfinal}]

Plot[Evaluate[vx[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

Plot[
      Evaluate[((ωb[
                    t] + λ (ϕ[t] - 
                        ArcTan[-x[t], -y[
                              t]])) /. λ -> 5) /. 
          sols], {t, 0, Tfinal}, PlotRange -> All]

Plot[Evaluate[ϕ[t] /. sols], {t, 0, Tfinal}, 
      PlotRange -> All]

gr5 = 
        ParametricPlot[Evaluate[{x[t], y[t]} /. sols], {t, 0, Tfinal}, 
          PlotPoints -> 1000, PlotRange -> All, 
          AxesLabel -> {"x", "y"}];

Show[gr1, gr2, gr3, gr4, gr5]